!=========Calcula período de um pêndulo por integral elíptica==========

PROGRAM exerB2

!======================================================================

IMPLICIT NONE

!========================Declaração de variáveis=======================

      INTEGER :: J
      REAL(8) :: f, h, mid, tetha0, integral, l
      REAL(8), PARAMETER :: g = 9.80665d0, delta = 0.0001

!===========================Leitura dos dados==========================

      WRITE(*,*) "Entre com o comprimento da haste l"
      READ(*,*) l
      WRITE(*,*) "Entre com o ângulo inicial tetha0"
      READ(*,*) tetha0
      WRITE(*,*) "Entre com o h para ser usado na quadratura"
      READ(*,*) h

!==================Faz a integral pela regra de Bode===================

      J = 0
      mid = -tetha0 + h! + delta
      integral = 0.0d0

      DO WHILE (mid <= (tetha0-h))
            !Calcula o meio do intervalo a ser usado como
            !como referencia no cálculo da área
            mid = -tetha0 + (2*J+1)*h! + delta
            !Regra de Simpson
            integral = integral + (h/3.0d0)*(f(mid+h, tetha0)+&
            &4.0d0*f(mid, tetha0)+f(mid-h, tetha0))
            J = J + 1
      END DO

!      DO
 !           mid = -tetha0 + J*4*h
 !           !Regra de Bode
 !           integral = integral + (2.0d0*h/45.0d0)*(7.0d0*&
 !           &f(mid, tetha0)+32.0d0*f(mid+h, tetha0)+12.0d0*&
 !           &f(mid+2.0d0*h, tetha0)+32.0d0*f(mid+3.0d0*h, tetha0)+&
 !           &7.0d0*f(mid+4.0d0*h, tetha0))
 !           J = J + 1
     !       write(*,*) integral
 !     END DO

!=======================Calcula o período e exibe======================

      WRITE(*,*) "O período do pêndulo é:", DSQRT(2.0d0*l/g)*integral

!======================================================================

END PROGRAM exerB2

!======================================================================

!==============================Funções=================================

REAL(8) FUNCTION f(tetha, tetha0)
      IMPLICIT NONE

      REAL(8) :: tetha, tetha0

      IF (DCOS(tetha) > DCOS(tetha0)) THEN
            f = 1.0d0/DSQRT(DCOS(tetha)-DCOS(tetha0))
      ELSE
            f = 0.0d0
      END IF
      write(*,*) tetha, f
END FUNCTION f

!======================================================================
